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Abstract 

v. 

An  investigation  of  the  spin-orbit  resonance  effects  about  the 
geosynchronous  orbit  was  undertaken  to  determine  the  existence  and 
feasibility  for  use  of  additional  stable  equilibrium  points  in  this 
regime  for  the  placement  of  communications  satellites.  The  Hamiltonian 
of  the  geopotential  was  developed  in  Delaunay  elements  using  first  and 
second  order  zonal  and  sectorial  harmonic  terms^  ( Jf  t  Ji »  J» » 

The  resulting  Hamiltonian,  which  was  valid  for  all  inclinations  and 
small  eccentricities,  was  reduced  to  a  single  degree  of  freedom  through 
a  series  of  transformations  to  allow  computer  generation  of  phase  por¬ 
traits  and  analysis  of  the  structure  and  librational  stability. 

Three  resonance  bands  were  discovered,  two  of  which  seemed  practi¬ 
cal  for  operational  use.  Although  the  second  of  these  bands  was  con¬ 
tained  within  the  primary  resonance  structure,  its  width  (246  meters) 
and  librational  period  (625  years)  appear  useful  for  satellite  placement 
at  the  additional  stable  equilibrium  points.  Additionally,  the  general 
technique  was  demonstrated  as  a  feasible  alternative  investigative  tool 

to  a  purely  analytic  approach.  . 
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I  Introduction 


With  the  growth  of  world-wide  communications  systems,  the  use  of 
geosynchronous  satellites  has  become  progressively  more  important  over 
the  past  decade.  Within  the  Department  of  Defense  the  geosynchronous 
communications  satellite  is  becoming  a  vital  link  in  C*  —  Command, 
Control,  and  Communications.  Unfortunately,  the  characteristics  of  the 
geosynchronous  orbit  present  two  major  difficulties. 

Because  current  communications  satellites  require  a  minimum  beam- 
width  separation,  the  fact  that  the  nominal  geosynchronous  orbit  is  a 
circular,  equatorial  orbit,  places  an  upper  limit  on  the  maximum  number 
of  satellites  which  can  inhabit  it,  resulting  in  it  being  considered  a 
scarce  resource.  Also,  when  basic  resonance  effects  are  considered, 
those  arising  from  the  triaxiality  of  the  earth's  figure,  it  is  discov¬ 
ered  that  this  orbit  has  only  four  primary  equilibrium  points,  with  only 
two  of  these  points  being  stable.  Therefore,  considerable  station¬ 
keeping  is  required  to  maintain  the  positioning  necessary  for  adequate 
communication  separation.  Modern  satellite  lifetimes  can  be  limited  by 
these  station-keeping  requirements. 

In  an  attempt  to  find  a  means  to  alleviate  these  problems,  an 
investigation  of  the  spin-orbit  resonance  effects  arising  from  the 
interaction  of  the  spin  of  the  nonspherical  earth  with  a  satellite  orbit 
was  conducted.  It  will  be  shown  that  these  resonance  effects  result  in 
multiple  equilibrium  points  at  near  geosynchronous  orbit  and  that  these 
equilibrium  points,  although  not  truly  geosynchronous,  may  be  useful  in 
permitting  an  expanded  use  of  this  orbit  for  communications  satellites. 


Background 


A  review  of  the  literature  reveals  that  an  extensive  amount  of 
research  has  been  done  concerning  the  effect  of  the  earth's  shape  on 
satellite  orbits.  The  majority  of  this  research  began  with  the  advent 
of  the  artificial  satellite  in  the  1950s.  In  particular,  the  effect  of 
the  earth's  oblateness  has  been  well  established  by  several  researchers. 
Many  varied  approaches  were  taken,  as  illustrated  in  the  November  1959 
issue  of  The  Astronomical  Journal  where  three  independent  treatments  of 
the  'main  problem'  of  satellite  theory  were  published. 

Selected  investigations  have  attempted  to  determine  the  effects  of 
the  earth's  longitude  dependent  harmonics  on  satellite  motion.  Because 
of  its  use  for  communications  satellites,  the  geosynchronous  orbit  has 
proven  to  be  of  great  interest. 

Blitzer  (Ref  7,8),  using  a  linearized  theory,  discussed  equilibrium 
solutions  for  a  geosynchronous  satellite  in  a  circular,  equatorial  orbit 
under  the  influence  of  the  principle  longitude  dependent  term,  Jzf 
Blitzer  (Ref  4,5)  later  treated  the  effect,  due  to  higher  order  tesseral 
harmonics,  on  a  geosynchronous  satellite  of  small  eccentricity  and 
inclination.  He  showed  that,  due  to  the  dominance  of  the  Ju  term,  only 
a  slight  displacement  of  the  equilibrium  positions  occurred,  although 
there  was  a  significant  change  in  the  librational  periods. 

Musen  and  Bailie  (Ref  22),  in  1962,  studied  the  motion  of  synchro¬ 
nous  satellites  incorporating  only  the  tesseral  harmonic  and  the  Jj 
and  zonal  harmonics.  Their  analytic  technique  agrees  with  Blitzer 
for  synchronous  motion  in  the  equatorial  plane. 

Allan  (Ref  1),  in  1965,  discussed  the  motion  of  nearly  circular  but 
inclined  synchronous  orbits.  In  1967  (Ref  2),  he  studied  the  effect  of 
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resonance  in  inclination  for  synchronous  satellites  in  nearly  circular 
orbits  when  the  positions  of  the  nodes  repeat  relative  to  the  rotating 
primary.  He  also  studied  the  effect  of  resonance  in  eccentricity  and 
inclination  for  a  synchronous  satellite  in  a  nearly  equatorial,  eccen¬ 
tric  orbit  which  occurs  when  the  longitude  of  the  line  of  apsides 
repeats  relative  to  the  primary  (Ref  3). 

In  a  book  published  in  1966,  Kaula  (Ref  19)  developed  expressions 
for  the  resonant  disturbing  function  of  the  geopotential  in  terms  of 
inclination  and  eccentricity  functions  and  derived  general  expressions 
for  the  variation  of  orbital  elements  due  to  arbitrary  zonal  or  tesseral 
harmonics . 

And,  in  a  series  of  articles  over  the  years,  Garfinkel  (Ref  15)  has 
developed  a  formulation  known  as  the  Ideal  Resonance  Problem  which  has 
found  wide  application  in  resonance  theory  by  treating  specific  cases  as 
perturbations  of  that  problem. 

From  the  preceding,  it  is  apparent  that  considerable  progress  has 
been  made  in  resonance  theory,  although  it  appears  at  this  time  that 
little  hope  exists  for  a  general  analytic  solution  except  for  certain 
specific  cases.  These  include  mostly  satellites  whose  mean  motion  is 
strictly  commensurate  with  the  earth's  rotation  rate  and  whose  orbits 
are  at  the  critical  inclination  or  have  zero  eccentricity.  It  should  be 
noted  that  while  Dallas  and  Diehl  (Ref  12)  claimed  such  a  general  solu¬ 
tion  in  1977,  they  were  later  shown  by  Jupp  (Ref  17)  to  have  made  a 
serious  error.  Garfinkel  (Ref  13),  in  his  1979  summary  of  the  Ideal 
Resonance  Problem,  lists  the  synchronous  satellite  with  non-zero  eccen¬ 
tricity  and  the  general  p:q  resonance  between  the  period  of  revolution 
of  the  satellite  and  the  rotation  of  the  earth  as  two  of  the  outstanding 


unsolved  problems  of  resonance  theory. 

As  a  result,  although  the  effort  continues  to  find  formal,  global 
solutions  to  resonance  problems,  most  recent  efforts  address  themselves 
to  specific  satellite  orbits  or  attack  the  problem  via  semi-analytic  or 
numerical  techniques  such  as  used  by  Nacosy  and  Diehl  (Ref  23). 

One  such  technique  is  the  computer  plotting  of  phase  portraits  of 
the  specific  system  Hamiltonian.  Although  this  method  has  been  used 
previously  to  provide  some  sort  of  physical  feel,  it  would  appear  that 
little  has  been  done  to  utilize  the  process  as  the  primary  investigative 
tool.  While  this  technique  provides  information  only  in  the  regime  of 
the  particular  inclination  studied,  it  is  valid  for  small  eccentricities 
and  permits  investigation  of  all  resonance  terms  associated  with  a 
particular  commensurability  ratio.  It  is  also  easily  modified  to  incor¬ 
porate  additional  harmonic  terms  or  study  different  commensurability 
ratios.  But  most  importantly,  it  circumvents  most  of  the  mathematically 
sophisticated  and  algebraically  laborious  techniques  required  in  the 
search  for  more  general  analytic  solutions. 


Objectives  and  Scope 

This  investigation  will  undertake  to  determine  whether  multiple 
equilibrium  points  exist  for  the  geosynchronous  resonance  and,  if  so, 
just  what  the  characteristics  of  these  equilibrium  points  are.  These 
characteristics  should  indicate  whether  further  investigations  are  jus¬ 
tified  to  determine  the  usefulness  of  these  points  for  the  placement  of 
communications  satellites. 


A  secondary  objective  of  this  investigation  is  to  demonstrate  the 


feasibility  of  the  general  technique  of  computer  generation  of  phase 
portraits  for  the  examination  of  geopotential  resonant  structures. 


The  scope  of  this  particular  investigation  will  be  limited  to  an 
examination  of  the  spin-orbit  resonances  arising  from  the  triaxiality  of 
the  earth's  figure.  The  contributions  from  the  first  and  second  order 
zonal  and  sectorial  harmonics,  the  J^,  (Jj)*,  J^,  and  terms  of  the 
geopotential,  will  be  examined  for  their  affect  on  a  near  geosjrnchronous 
orbit. 


General  Approach 

The  Hamiltonian  of  the  geopotential  will  be  developed,  using 
Delaunay's  elements,  by  combining  the  secular  first  and  second  order 
zonal  harmonic  terms  with  the  nominal  two-body  Hamiltonian,  as  extracted 
from  Brouwer's  1959  article  in  The  Astronomical  Journal  (Ref  9),  and 
then  developing  the  perturbing  function  arising  from  the  primary  sector¬ 
ial  harmonic  term,  J,!.  Once  obtained,  this  Hamiltonian  will  then  be 
converted,  by  means  of  a  contact  transformation,  to  another  set  of 
elements  which  will  facilitate  work  for  the  geosynchronous  case.  At 
this  point,  each  resonance  term  will  be  treated  independently  and 
another  canonical  transformation  will  be  used  to  allow  reduction  to  a 
single  degree  of  freedom.  This  form  will  allow  the  resonance  effects  to 
be  analyzed  by  use  of  phase  portraits  of  this  single  degree  of  freedom. 
From  these  phase  portraits,  stable  equilibrium  points  can  be  determined 
and  the  stability  analyzed. 


Sequence  of  Presentation 

Chapter  II  presents  the  theoretical  development  of  the  system 
Hamiltonian,  the  equations  of  motion  for  each  resonance  term,  and  the 
methods  for  determining  librational  stability.  Chapter  III  discusses 
the  computer  application  of  these  theoretical  developments  and  Chapter 
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IV  presents  an  analysis  of  the  results  and  conclusions  of  the  study  as 
well  as  recommendations  for  further  research. 


II  Theory 


Geopotential  Hamiltonian 

In  determining  the  spin-orbit  resonance  effects  for  a  near  geosyn¬ 
chronous  orbit,  initially  the  Hamiltonian  for  the  geopotential  must  be 
obtained.  The  general  form  for  the  Hamiltonian  may  be  written 

F  -  F.  t  F.  »  P._  .  F.,  »  R., 

where  F,  represents  the  contribution  of  the  unperturbed  potential,  Fj, 
F,^,  and  F«^  are  the  secular  contributions  of  the  Jj ,  Jf*  and  terms, 
and  R,2  is  the  disturbing  function  resulting  from  the  J,,  term.  The 
first  four  terms  may  be  extracted  from  a  1959  Astronomical  Journal  arti¬ 
cle  by  Brouwer  (Ref  9)  and  written  ignoring  long  period  variations  as 


1  3  H 

-  -p^(- 1  + 


1 +  (-)'*) 
4L‘*  ‘32'G'^  ^  S'G"^  'g'  ' 
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where 


L  ■  /pa 
G  -  L/(l-e*) 
H  ■  G(cosi) 


£  >  mean  anomaly 

g  ~  argument  of  perigee 

h  »  longitude  of  the  ascending  node 


and  u  is  the  geocentric  gravitational  parameter,  is  the  geocentric 
mean  equatorial  radius,  a  is  the  semi-major  axis,  e.  is  the  eccentricity, 
and  i  is  the  inclination. 

The  general  form  of  the  geopotential,  as  listed  in  Hagihara  (Ref 
16:459),  is  now  used  to  determine  the  disturbing  function  due  to  the 
principal  sectorial  harmonic,  J^t* 

hR* 

Rj2  =  -^J**P|(sin8)sin(2X-21t  j  ) 

Here  B  is  the  latitude,  X  is  the  longitude,  and  X^j  is  the  longitude 
associated  with  the  coefficient.  Since 

P|(sin0)  “  3cos*B 
l»R* 

R*2  *  -^J**(3cos*B)sin(2X-2Xj, ) 


Figure  1.  Orbit-Bquator-Meridian  Triangle 

To  express  the  disturbing  function  in  terms  of  the  Delaunay  ele¬ 
ments,  B  and  X  must  first  be  converted.  Examination  of  Figure  1  indi¬ 
cates  that  the  right  ascension,  a,  of  a  satellite  at  point  P  would  be 
expressed  as 


o  ■  h  +  B 

and  since  the  satellite's  longitude  is  the  difference  between  its  right 


ascension  and  the  Greenwich  sidereal  time,  6,  then 

X»a-e“h  +  6-  e 

The  Greenwich  sidereal  time  may  also  be  expressed  as 

0  =  not 
so 

X  =  h  +  «  -  n,t 

Letting 

if  =  n,t  +  Xjo 

so  as  to  simplify  the  following  expressions,  and  performing  some  trigo 
nometric  expansions 

3cos*B8in(2h-2^+26)  =  3cos*0  [sin(2h-24>)cos26  +  cos(2h-24))sin26] 

cos26  =  cos*4  -  sin*8 
sin26  2sin8cos8 
From  spherical  trigonometry, 

sinfi  *  sinucosi/cos8 
cos8  =  co8u/cos0 

where  u  is  the  argument  of  latitude,  so 

cos26  =  (cos*u  -  sin*ucos*i)/cos*0 
sin28  *  2sinucosucosi/cos*0 

and  therefore 

3cos*P8in(2h-2t+26)  ”  38in(2h-2t)(co8*u  -  8in*uco8*i) 

+  6co8(2h-2ii)(8inucosuco8i) 

Note  that  the  latitude  dependence  of  the  disturbing  function  has  now 
disappeared  and  it  is  no  longer  necessary  to  express  6  in  terms  of  the 
Delaunay  elements. 


Now,  by  substituting 


1  1^1 
co8*u  '2  ■2^0820 


•  »  11, 

8in*u  "2  “  -^00820 


1  •  o 

8inuco8u  ®  ■28in2u 


into  the  above  equation  and  amplifying 


3co8  *  3  a  in  (  2h-2  tJi+2  5  ) 


1+C08  *  i  )  8  in(  2h-2  )  cos  2u 
+  ^in*  isin(2h-2t|i) 

+  3cosicos(2h'*2i)sin2u 


Further,  noting  that 


8in(2h--2i|i)co82u  =  -|■[8in(2h-2♦+2u)  +  sin(2h-21*-2u)] 
co8(2h-2<(i)8in2u  =  •j(8in(2h-2^+2u)  -  8in(2h-2t|»-2u) ] 


and  substituting  these  expressions  into  the  last  equation  yields 


3co8*68in(2h-2t|i+26) 


^(l+co8i)*8in(2h-2^+2u) 


+  -^sin*  i8in(2h-2+) 


+  ^(l-co8i)*8in(2h-2^-2u) 


Finally,  writing 


u  «  f  +  g 

where  f  is  the  true  anomaly,  and  substituting  back  into  the  expression 
for  the  disturbing  function  results  in 

“^e  3 

R,,  -  -^J,,[^(l+co8i)*8in(2g+2h+2f-2*) 
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+  ■^in*isin(2h-2t) 

+  ^(l-co8i)*sin(-2g+2h-2f-2t)] 

To  express  this  term  solely  as  a  function  of  the  Delaunay  elements,  r 
and  f  must  be  expanded  in  terms  of  the  mean  anomaly,  1,  and  substituted 
into  the  above  expression.  Note  that  up  to  this  point  in  the  develop¬ 
ment  of  the  geopotential  Hamiltonian,  a  completely  general  expression 
has  been  formed,  valid  for  all  eccentricities  and  inclinations  and 
limited  only  by  the  number  of  zonal  and  sectorial  harmonic  terms 
included  in  the  treatment. 

Using  r  and  f  as  listed  by  Brouwer  and  Clemence  (Ref  10:76-77)  and 
retaining  terms  to  the  seventh  power  in  eccentricity, 

—  ®  1  +  e(-  cost)  +  ~  •^0821)  +  e*(4co8l  -  TCOsSt) 

a  z  2  o  o 

+  e*(-|co82l  -  ^o84a)  +  e*(-  ~ 

1  2  27 

+  e*(-  y^082i  +  ycosAt  -  -^oaSl) 

.  1  567  .  4375  ..  16807 


f  =  t 


+  e'C-  •^8in2)l  +  -^^8104)1)  +  e’C-sl^sinl  -  T^inSl  + 
24  96  96  64 

^  17  .  451  ^  1223  .  ,.v 

+  e‘(y^sin2i  -  ^8in4A  +  -^in6l) 

,,  107  ..  .  95  .  5957  ....  47273  . 

^  ®  *  5Y281"3A  -  ^8in5i  +  32^1^71) 


Now,  converting  to  canonical  units  (R^  ■■  p  ■■  1)  and  substituting  these 


expressions  into  R. . 


R*»  ■ 


*^**3  11  S  lA'i 

-jr<;^(l+coai)*I(-  ^  “  384«*  ••  ^8432®^^®^“^ 

+  (-  ^ge  -  yggfi*  "  30720® ^  ^ A-2g-2h+2 ♦) 

+  (1  -  |e*  +  ||e'  -  ~|e*)8in(21+2g+2h-2*) 

^  ^  fir*  -  ^^)^i-(3U2g*2l.-2^) 

81  81 

^  1280®’  ■  20^'>«i“(34-2g-2h+2*) 

+  (iZe*  -  ig^e-  +  ^||e*)8in(44+2g+2h-2<») 


+  (-  ^*)8in(4l-2g-2h+2*) 

*  -  W^'  *  ^|f|feOsin(5t.2g.2h-2*) 


*  "  ^^Osin(5t.2g.2h-2 

*  iMfif®  ^  ^ g“2h»2 » ) 

+  (-^e^  -  i|||4*)8in(6l+2g+2h-2t) 

^  ,228347  ,  3071075,^  .  ^  . 

^  3840  ®  18432  ®  )8in(7t+2g+2h-2<i) 

+  (- y  l^-e  * )  8  in  ( 8it-»-2g»2h-2  ») 

+  (--y||3^^e08in(9i->-2g42h-2»)] 

+  |8in*i[(l  +  |e*  +  +  ||e*)8in(2h-2*) 

^  ^  fl®’  ^  ffl®‘  ^  l|f|f4’)8in(U2h-2*) 


*  (.  la  .  12L*  IIL,*  1^309^ 

'  2  12^  6144 


4  (|e»  +  |e*  + 

4  4  64 


W  *  ^*)8in(2l42h-2*) 
4  64 


~||e')8in(l-2h42*) 
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+  (^‘)sin(4i-2g+2h-2i|») 

.  ,  845  ,  .  32525  ,  208225  n 

48®  768  ®  ~  6144  ®  )®i“^5«.+2g-2h+2+) 

^  (Y|fff|e')«i«(51-2g+2h-2*) 

+  ^^^e*)8in(6A+2g-2h+2i|.) 

lo  iou 

3840  ®  18432  ®  2g-2h  2(|i) 

+  (-  ^|||^*)sin(8A+2g-2h+2^) 


+  (-  ^ ^ n tsO^® ^  ) 8 in(  9 t^-2g-2h-^2 » )  ]  > 


Although  this  form  of  8,2  now  limited  to  small  eccentricities,  the 
practical  limit,  as  noted  by  Kovalevsky  (Ref  20:55),  is  of  sufficient 
magnitude  as  to  have  no  further  affect  on  this  treatment.  Obviously, 
for  small  eccentricities,  carrying  the  expansion  to  the  seventh  power  in 
eccentricity  should  allow  sufficient  accuracy  for  this  investigation. 

Examination  of  the  above  expansion  of  Rjx  indicates  the  existence 
of  many  resonance  terms,  permitting  investigation  of  resonances  other 
than  the  geosynchronous  resonance  merely  by  selecting  the  appropriate 
terms.  For  the  geosynchronous  case,  however,  resonance  occurs  when  the 
mean  motion  of  the  satellite  is  commensurate  with  the  mean  rotation  rate 
of  the  earth,  leading  to  selection  of  terms  including  2t-2i(i  in  the 
trigonometric  argument. 

Retaining  only  the  near  geosynchronous  resonance  terms  and  comple¬ 
ting  the  transformation  to  Delaunay  elements  yields 


X 


^2  2  -J 


-)*(- 
g'  ' 


233 

288 


96^L^ 


+  ■j||(^)*)sin(2A+2g+2h-2*) 


+ 


+  31.G., 

G  240  240  X  240^l' 


■^(|)  *  )8in(21-2g+2h-2t)  ] 


The  geopotential  Hamiltonian  is  now  expreaaed  completely  aa  a  function 

of  the  Delaunay  elements  and  time  as 

F  =  F(l,g,h,L,G,H,t)  =  F,  +  F,  +  F,  +  F.  +  Rj, 

8  8 


Canonical  Transformation  for  the  Geosynchronous  Case 

In  order  to  simplify  calculations  later  in  this  development,  a  new 
set  of  variables  is  selected  which  comprises  a  canonical  transformation 
such  that  two  of  the  new  momenta  variables  are  small  for  the  geosynchro¬ 
nous  equatorial  case.  These  new  momenta  are  defined  as 

X  *  L  =  ^ 

Y  =  L  -  G  =  /vail  -  /(!l-e»)] 

Z  =  G  -  H  =  /va( l-e*)(l-cosi) 

As  a  result,  the  generating  function  becomes 

Si  =  X£  -►  (X  -  Y)g  (X  -  Y  -  Z)h 
so 


x  =  i  +  g-^h  y  =  -g-h  z  =  -h 
The  Hamiltonian  for  the  geopotential  may  now  be  expressed  as 
F  =  F(x,y,z,X,Y,Z,t)  =  F,  i-  F,  F,^  +  F,^  R*, 

where 


F. 


1 

2X* 


J* 


+  fB*) 


F, 


+  |a‘(4  -  24B  +  48B*  -  36B»  +  9BM 

o 

-  -IfA’C-  8  +  32B  -  44B*  +  28B*  -  78")! 


R*a  =  ^[|(4-4B+B*)(1  -  5C  +  ^C*  - 

"  ^  ^‘)sin(2x-2*) 

+  |(2B-B»)(|c  + 

-  ’  ^•)8in(2x+2y-2z-2*) 

*  !“’<6<='  * *  T^’  -  ■S5C*).in(2x*4,-4z-2*)] 


Note  that  for  the  near  geosynchronous  case,  the  coefficients  B  and  C  are 
small  while  the  coefficient  A  is  approximately  unity.  From  this  obser¬ 
vation  it  can  be  seen  that  the  first  resonance  term  is  the  strongest, 
followed  by  the  second  and  third  terms  respectively. 


Resonance  Transformations 

Now  that  a  general  form  of  the  geopotential  Hamiltonian  has  been 
developed  for  the  geosynchronous  case,  each  of  the  three  resonance  terms 
may  be  individually  examined.  That  is,  the  geopotential  Hamiltonian  may 
be  considered  to  consist  of  the  first  and  second  order  zonal  terms  along 
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with  a  single  resonance  term  from  the  disturbing  function.  The  justifi¬ 
cation  for  this  approach  lies  in  the  assumption  that  in  the  immediate 
region  of  each  resonance  band,  the  equations  of  motion  are  dominated  by 
the  corresponding  resonance  term.  The  relative  strengths  of  the  three 
terms  should  guarantee  this  result. 

To  facilitate  this  examination,  phase  portraits  of  the  geopotential 
Hamiltonian  will  be  used.  In  a  phase  portrait,  each  generalized  coordi¬ 
nate  and  its  associated  momentum  make  up  two  of  the  dimensions  of  the 
phase  space.  The  real  advantage  to  this  representation  is  that  paths 
corresponding  to  particular  unique  solutions  of  the  system  Hamiltonian 
will  be  produced  (Ref  21:173). 

At  this  point,  however,  the  geopotential  Hamiltonian  has  six  dimen¬ 
sions  arising  from  the  three  generalized  coordinates  and  their  conjugate 
momenta.  And  since  time  is  still  explicit  in  the  Hamiltonian,  an  addi¬ 
tional  dimension  is  introduced.  These  seven  dimensions  constitute  what 
is  referred  to  as  a  motion  space.  Because  it  is  extremely  difficult  to 
conceptualize  a  trajectory  in  a  seven-dimensional  space,  it  would  be 
extremely  advantageous  if  the  geopotential  Hamiltonian  could  be  limited 
to  only  two  variable  dimensions  to  permit  plotting  of  these  trajecto¬ 
ries.  This  result  may  be  accomplished  by  reducing  the  geopotential 
Hamiltonian  to  a  single  degree  of  freedom  and  eliminating  the  explicit 
time  dependence  through  the  use  of  a  canonical  transformation. 

In  each  case,  a  new  canonical  transformation  will  be  defined  such 

that 

*  * 

F  -  F  (_,_,s,Q,R,S,_)  -  F(x,y,z,X,Y,Z,t)  -  — 
where  the  associated  generating  function  will  be  of  the  form 


Sj  *  S{ (x tY , z , Q,R,S I t ) 

These  transformations  will  yield  a  Hamiltonian  for  each  resonance  of  the 
form 

F*  =  F,  +  F,  +  F,^  +  F»^  +  Rja  +  n,S 

wi  th 


•^*13  3 

^  -  |b  >  |b*) 


's  4S‘*^32 


rr[^A*(-  I  +  ^B  +  -^B*  -  4B*  +  B') 


+  |a*(4  -  24B  +  48B*  -  3tB*  +  9B*) 


15,, 


32 


■A’(-  8  +  32B  -  44B*  +  28B*  -  78“)] 


F, 


These  transformations  will  be  done  in  such  a  way  that  the  functional 
forms  of  Fj,  F,,  F,^,  and  F,^  will  always  be  as  listed  above.  The  only 
change  from  resonance  to  resonance  will  be  the  functional  dependence  of 
the  coefficients  A,  B,  and  C  on  the  generalized  momenta  Q,  R,  and  S. 

The  particular  resonance  term  selected  from  the  disturbing  function  will 
be  transformed  so  that  its  trigonometric  argument  will  depend  upon  a 
single  generalized  coordinate. 

By  selecting  a  transformation  of  this  form,  three  constants  of  the 
motion  result.  Since  neither  q  or  r  appear  explicitly  in  the  resulting 
geopotential  Hamiltonian,  Q  and  R  are  constants.  Also,  because  the 
geopotential  Hamiltonian  describes  a  natural  system,  it  represents  the 
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total  system  energy,  which,  due  to  the  lack  of  an  explicit  time  depend¬ 
ence,  is  also  conserved.  As  a  result,  the  new  Hamiltonian  is  reduced  to 
one  degree  of  freedom,  permitting  the  use  of  phase  portraits. 


For  the  first  and  primary  resonance  condition, 

R»  =  ;^(4-4B+B*)(l  -  5C  +  -  Ifc*  ^*)sin(2s) 

By  noting  that  the  trigonometric  argument  of  this  term  expressed  as  a 
function  of  x,  y,  z,  and  t  (as  developed  in  the  previous  section)  was 

2x  —  2nQ t  —  2X22 


s  ™  X  —  n®  t  —  X  2  2 

By  letting  q  =  y  and  r  =  z,  the  generating  function  becomes 


and  therefore 


Sa  *  Qy  Rz  S(x  -  net  -  X22) 


X*S  Y  =  Q  Z=R 


A  =  B  =  C  *  ^ 

S-Q  S-Q  S 

Equilibrium  points  exist  where  the  time  derivatives  of  the  only 
remaining  generalized  coordinate  in  the  Hamiltonian  and  its  associated 
momentum  are  both  zero.  That  is  when 

ds  3F*  -  .  dS  3F*  ^ 

dt  3S  dt  3s 

Therefore,  for  the  first  resonance  term,  the  equilibrium  points 
exist  where 


19 


_  *  rA»(15  .15  «  ,  ^  75  ,  _  75 

1  ^  2®  8®  8®  32®  ^ 

+  A*(-  I  +  27B  -  +  69B*  -  -^B') 

.  A»(-  ^  .  108B  -  l^B*  .  i2|73.  _  13|53.^ 

+  A*(-  +  120B  -  +  -^B*  -  -m^B"*)] 

-  ^[(A*-AM(-  ^  +  75B  -  ^B*  .  ^B*  -  ^BM 

+  A*(-  ^  +  90B  -  ^B*  +  210B*  -  ^B**) 

+  A*(^  -  200B  +  ^^B*  -  ^B*  +  -m^B")] 

2.  4  2  lo 

-  ^^[(4-4B+B*)(-  6  +  35C  -  46C*  +  38C*  -  ^Hc**  +  -^Hc*  -  ||c*) 

+  A(4B-2B*)(1  -  5C  +  -  ||c*  +  ^*)]8in(28) 

-  n,  =  0 


f  -  ^(4-4B.B*)(1  -  5C  .  -  Ifc*  .  1||c^  -  ||c*  .  ^*)co8(28) 


For  the  8econd  re8onance  term. 


*  i|c- *  2|c- 


395^>  ^  423^,  141..,  .  , 

*  ~3^  ■  -6^*)8in(28) 


and  chooaing 


8»x+y-2-  n,t  -  X,j 


q  ■  y  r  “  a 


yielda 


§2  “  Qy  +  Bz  +  S(x  +  y  -  2  -  n,t  -  X,,) 
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and  therefore 


X=S  Y  =  Q  +  S 


A  =  -  1 

^  Q 


Now,  the  equilibrium  condition  becomes 


ds  1 


fA*(-  I  +  |b  -  |b*)  +  A*(4  -  h 


+  A*(15 


-  2sr-((2B-B*)(^  -  22C  -  ^ 

-  A(2-2B)(|c  +  -  ■^*)]8in(28) 

-  n,  »  0 


g  -  ^(2B-B*)(|c  +  i|c*  .  -h  -  1|^.)co8(28)  -  0 


For  the  final  resonance  term. 


-  "  T!^*  -  "  4^’ 


2lo^‘^8in(28) 


8  -  X  +  2y  -  2z  -  n,t  -  X„  q  -  y  r  » 
S.  ■  Qy  ♦  R*  +  S(x  +  2y  -  2z  -  n.t  -  X,,) 


X  =  s 


Y  -  Q  +  2S 


Z  “  R  -  2S 


A  =  - 


2S+Q 

S 


a 


and  the  resonance  condition  becomes 


-  5T  -  >•(-  I  ♦  i.  -  !»■)  ♦  A<(f  -  9B  . 


-  15,  .  45,,  .  75  .  _  75 
4s7t1*  1-4  -  -2“  8*  *  8®  32'  ’ 

.  *.(-  5i  ^  8i,  _  ,  m,.  _  2i7,., 

»  A’(“  -  90B  ♦  ^^B‘  -  ♦  -^A*) 


»  -  ^B  *  ili5„.  -  3U,.  ,  1^B>)1 

^KA.-A-K-  if  .  75B  -  5^8.  *  ^B>  - 


+  A*(~  -  225B  +  +  ^B") 

2  4  2  16 

+  A*(-  +  425B  -  ^^B*  +  525B*  -  -^tPb')] 

2  4  16 

-  1^.  -  lie.  *  l|c*  -  |1C.  . 


-  A(4B-2B.)(ic.  * 

^  n,  *  0 

ft  ■  *  T5<=‘  -  #=•  *  4^'  -  So®’'"-”*’  ■  “ 

Obviously,  for  all  three  resonances,  the  equilibrium  points  will  occur 
where 


(2n+l  )ir 


n  •  0,1 ,2 ,3 


This  observation  can  be  made  from  a  brief  examination  of  the  first  and 
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second  partial  derivatives  of  the  geopotential  Hamiltonians.  The  gen¬ 
eral  form  of  the  first  partial  derivatives  is 

* 

- —  “  kco8(2s) 

as 

where  k  is  a  constant  for  some  fixed  value  of  S  corresponding  to  the 
value  at  the  equilibrium  point.  Now,  the  second  partial  derivatives 
have  the  form 

7-  *  -  2ksin(28) 

08 

indicating  a  negative  value  when  n  is  even  and  a  positive  value  when  n 
is  odd.  Because  a^F  /as*  is  of  order  1/S^,  which  is  always  greater  than 
zero,  then,  by  the  Second  Derivative  Test,  a  local  minimum  exis  when  n 
is  odd  and  a  saddle  point  exists  when  n  is  even.  The  local  minima 
correspond  to  stable  equilibrium  points  and  the  saddle  points  to  unsta¬ 
ble  equilibrium  points.  This  treatment  assumes  that  the  cross  partial 
terms  are  negligible,  which,  due  to  the  system  geometry,  should  be 
expected.  This  assumption  was  subsequently  verified  through  the  use  of 
numerical  techniques. 

Although  each  resonance  term  has  identical  values  of  s  defining  the 
equilibrium  points,  their  interpretations  are  slightly  different.  For 
the  first  term, 

s  =  I  +  g  +  h  -  n,t  -  X,, 

which  for  the  nearly  circular,  nearly  equatorial  case  of  interest  to 
this  investigation 

ft  s  f 
ft  +  g  z  6 

which  means  that  the  stable  equilibrium  points  exist  where  the  'mean 
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longitude'  of  the  satellite,  X  ,  is 

n 


,  (2n+l)»  .  . 

V  •  —4— 


n  =  odd 


Since  the  values  of  s  for  the  other  two  resonance  terms  vary  only 
by  subtracting  differing  multiples  of  the  argument  of  perigee,  g,  the 
stable  points  for  these  resonance  bands  exist  where 


,  (2n+l)x  .  .  , 

•  — * — 


n  =  odd  (Second  Resonance) 


^  jj  (2n>l )t  +  2g  +  Xj,  n  =  odd  (Third  Resonance) 

m  4 

Unstable  points  occur  for  even  values  of  n. 

The  corresponding  values  for  S  are  somewhat  more  difficult  to 
determine.  Computer  analysis  will  be  used  to  determine  these  values  and 
to  plot  the  structure  of  the  resonances. 

Librational  Analysis 

Once  the  stable  equilibrium  points  have  been  located  for  each 
resonance,  there  are  basically  two  methods  for  obtaining  values  for  the 
period  of  libration  about  them.  The  first  method  would  involve  the 
evaluation  of  a  closed  line  integral  along  a  contour  of  constant  energy. 
Since  the  equations  of  motion  were  developed  for  the  remaining  general¬ 
ized  coordinate  and  momentum,  the  period  of  libration  for  any  given 
energy  contour  may  be  expressed  as  the  closed  line  integral  of  either 

(-J  Ot  /- 


where  §  and  s  are  the  time  derivatives  of  S  and  s,  respectively,  and  are 
each  functions  only  of  S  and  s.  To  attempt  to  determine  the  period  of 
libration  analytically  using  this  method  would  be  extremely  difficult. 


and  yet,  a  numerical  evaluation,  while  computationally  intensive,  should 
yield  very  good  results. 

The  second  method  involves  the  use  of  linear  analysis  to  estimate 
the  basic  oscillatory  period  in  the  immediate  region  of  the  stable 
equilibrium  points.  While  not  yielding  as  specific  a  result  as  the 
previous  method,  it  is  much  less  computationally  intensive  and  should 
result  in  an  answer  sufficient  for  the  purposes  of  this  investigation. 

Using  the  Chain  Rule  of  calculus 
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Now,  since 
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therefore 
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Assuming  again  that  the  cross  partial  terms  are  small  compared  to  the 
other  terms,  the  resulting  system  of  equations  becomes 
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or 


58  =  -  ii),5S 
5§  =  tii^Ss 

(uj  and  a>2  are  constants,  as  all  partial  derivatives  are  evaluated  at 
the  equilibrium  points).  Taking  the  time  derivative  of  the  first  equa¬ 
tion  and  substituting  into  it  the  second  equation  yields 

58  =  -  ((i)|Ii)2)58 

which  is  the  equation  of  motion  for  a  harmonic  oscillator  with  frequency 


/uiUj . 

While  this  method  is  somewhat  more  involved  analytically,  the 
resulting  frequency  should  prove  very  easy  to  evaluate  numerically. 
Using  an  approximation  of  the  limit  definition  of  a  second  partial 
derivative  and  the  fact  that  the  first  partials  are  identically  zero  at 
the  equilibrium  points  yields 


u.  a 


3F  (S+AS.8)/8S 
AS 


3F  (S.8-t-As)/3S 
As 


Since  these  partials  have  already  been  calculated  above,  it  is  a  simple 
matter  (once  these  equations  are  programmed)  to  approximate  the  fre¬ 
quency  of  oscillation  by  using  small  values  for  As  and  AS. 
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Ill  Computer  Application 


In  implementing  the  theory  described  in  Chapter  II,  four  basic 
programs  were  developed.  This  chapter  will  describe  those  computer 
algorithms  used  for  the  determination  of  resonance  values,  plotting  of 
the  phase  portrait,  and  the  two  methods  of  determination  of  the  libra- 
tion  period  about  the  stable  equilibrium  points.  Listings  of  the 
computer  algorithms  employed  are  provided  in  Appendix  A. 

Crucial  to  the  implementation  of  all  the  programs  used  in  this 
investigation  was  the  programming  of  the  equations  of  motion  developed 
in  Chapter  II.  These  equations  were  located  in  the  subroutine  FDF  and, 
when  passed  the  values  of  s  and  S  along  with  the  value  of  the  resonance 
of  interest,  provided  the  values  of  the  Hamiltonian  and  its  two  partial 
derivatives. 

Resonance  Value  Determination 

The  initial  step  taken  to  ascertain  the  utility  of  each  resonance 
structure  for  use  in  the  placement  of  satellites  was  to  find  the  loca¬ 
tion  of  the  equilibrium  points,  both  stable  and  unstable,  and  the  width 
of  the  resonance  structure.  The  program  VALUE  was  written  for  this 
purpose. 

Since  the  equilibrium  values  of  s  were  readily  determined  from  the 
equations  of  motion,  finding  the  location  of  the  equilibrium  points 
required  only  that  the  value  of  S,  for  the  zero  point  of  the  time 
derivative  of  s,  be  determined  along  these  equilibrium  values  of  s.  To 
accomplish  this  goal,  a  binary  search  for  the  zero  point  of  the  time 
derivative  was  employed  in  the  subroutine  RESON.  Once  the  equilibrium 
values  for  S  had  been  determined,  the  location  of  the  four  equilibrium 


points  for  each  resonance  structure  would  be  known.  To  give  more  of  a 
physical  feel  for  their  location,  their  radial  distance  from  the  nominal 
two-body  geos3mchronous  orbit  was  also  determined. 


a 


The  width  of  the  resonance  structure  was  found  by  determining  the 
roots  of  the  Hamiltonian  at  the  value  of  the  unstable  equilibrium  points 
along  the  azimuths  of  the  stable  equilibrium  values  of  s.  Initially, 
the  function  values  of  the  geopotential  Hamiltonian  for  the  stable  and 
unstable  equilibrium  points  were  determined.  Since  the  function  value 
of  the  unstable  equilibrium  points  defines  the  boundary  of  the  stable 
region,  the  roots  of  these  values  along  the  azimuth  of  the  stable  equi¬ 
librium  points  describes  the  maximum  width  of  the  resonance.  In  the 
subroutine  ROOTS,  Newton's  method 

*i+l  “  *i  ' 

for  finding  the  roots  of  an  equation  was  employed.  Once  the  two  roots 
were  determined,  the  distance  between  them  was  found  by  converting  back 
to  distance  units  (DUs)  and  finally  to  meters. 

The  results  of  VALUE  should  serve  to  give  an  idea  of  the  physical 
location  of  the  equilibrium  points  along  with  a  feeling  as  to  whether 
the  size  of  the  resonance  structure  is  adequate  for  satellite  placement. 


Plotting  of  Phase  Portraits 

Once  the  general  characteristics  of  a  given  resonance  structure  are 
determined,  the  actual  structure  of  the  resonance  may  be  plotted  as  a 
phase  portrait.  It  should  first  be  noted,  that  since  the  basic  struc¬ 
ture  of  each  of  the  three  resonance  bands  considered  in  this  investiga¬ 
tion  is  identical,  the  technique  used,  with  some  modifications,  is  the 
same  for  all  three  bands. 

In  order  to  generate  the  phase  portraits,  a  contour  following 
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subroutine,  furnished  by  Dr.  William  E.  Wiesel  and  modified  slightly  by 
the  author,  was  used.  This  subroutine,  CONTUR,  which  utilizes  Newton's 
method  for  finding  roots  generalized  to  functions  with  two  independent 
variables,  was  employed  in  conjunction  with  a  plotting  routine  for  the 
CALCOMP  1038  plotter,  in  the  subroutine  DRAW,  in  generating  the  phase 
portrait  in  the  program  RSPLOT. 

Librational  Periods 

Each  of  the  two  methods  described  in  Chapter  II  for  determining  the 
period  of  libration  about  the  stable  equilibrium  points  was  implemented. 
In  the  application  of  the  line  integral  approach  in  the  program  TIME, 
the  subroutine  DRAW  was  modified  to  create  the  subroutine  PERIOD.  This 
approach  was  used  in  order  to  take  advantage  of  the  already  developed 
algorithm  for  generating  the  points  for  a  given  energy  contour.  As 
these  points  were  generated,  an  elementary  numerical  integration  tech¬ 
nique  was  employed  to  determine  the  value  of  the  integral.  To  simplify 
the  calculation  of  the  distance  between  successive  points,  an  approxima¬ 
tion  was  made  to  use  the  largest  of  the  two  step  sizes  of  s  or  S.  This 
method  is  actually  an  adaptation  of  the  technique  used  in  CONTUR  to  find 
the  next  point  along  the  contour.  Switching  of  which  step  size  was  used 
necessarily  meant  switching  of  which  form  of  line  integral  was  being 
used,  but  this  switching  was  also  easily  performed  using  the  existing 
structure  of  CONTUR.  Once  the  value  of  the  line  integral  was  deter¬ 
mined,  the  result,  expressed  in  time  units  (TUs)  was  converted  to  mean 
solar  days. 

The  second  approach,  using  the  method  of  linearization,  proved  as 
easy  to  implement  as  expected.  In  the  program  TIMES,  values  of  the 
second  partials  of  the  Hamiltonian  with  respect  to  both  s  and  S 


(3*F  /38*  and  3*F  /3S*)  were  calculated  numerically  from  the  values  of 
the  first  partials  provided  through  the  subroutine  FDF.  These  values 
were  then  multiplied  together  and  the  square  root  taken  to  provide  the 
frequency  of  oscillation.  This  frequency  was  then  converted  to  a  period 
in  mean  solar  days. 
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IV  Analysis 


With  the  design  of  the  computer  algorithms  completed,  the  remaining 
portion  of  this  investigation  now  centers  on  an  analysis  of  methods  used 
and  their  validity  and  implications.  Before  turning  to  the  actual 
examination  of  the  results,  it  is  first  necessary  to  stipulate  the 
assumptions  made  in  their  generation. 

AsstMuptions 

Although  the  theory,  as  developed  in  Chapter  II,  for  the  investiga¬ 
tion  of  the  geosynchronous  spin-orbit  resonances  is  quite  general, 
certain  assumptions  had  to  be  made  during  the  computer  application  of 
this  theory  in  order  to  obtain  the  results.  None  of  these  assumptions, 
however,  will  cause  any  significant  changes  in  the  conclusions  to  be 
drawn  from  this  study.  Throughout  the  development  of  the  computer 
algorithms,  an  effort  was  made  to  permit  these  assumptions  to  easily  and 
readily  changed  to  allow  further  investigation  of  this  subject. 

Necessarily,  a  certain  set  of  constants  to  be  used  for  definition 
of  the  geopotential  was  required.  These  constants,  as  derived  from  the 
SAG  Standard  Earth  III  (Ref  13),  are  listed  in  Table  I  below. 

TABLE  I 

SAG  Standard  Earth  III  Constants 

J,  -  1082.637  X  10"* 

J»  -  -  1.617999  X  10"* 

J,,  -  2.7438636  x  10"* 

n,  -  5.86729371  x  10"*  rad/TU 

R  » 
e 


6.378140  X  10*  meters 


Another  assumption  involved  the  setting  of  the  values  of  the  two 
constant  momenta,  Q  and  R.  In  the  development  of  the  computer  algo¬ 
rithms,  it  was  assumed  that  it  would  be  advantageous  to  initially 
specify  some  particular  nominal  eccentricity  and  inclination  for  inves¬ 
tigation.  Careful  examination  of  the  equations  of  motion  developed  in 
Chapter  II  along  with  their  implementation  for  plotting  of  the  phase 
portrait  in  Chapter  III  indicated  that,  as  the  satellite  moves  along  a 
contour  of  constant  energy,  the  requirement  that  Q  and  R  are  constant 
dictates  that  the  eccentricity  and  inclination  be  constantly  changing. 
Since  these  changes  should  be  small,  the  values  of  Q  and  R  were  deter¬ 
mined  by  using  the  nominal  values  of  eccentricity  and  inclination  of 
interest  at  the  nominal  two-body  geosynchronous  radius.  Admittedly, 
this  determination  is  somewhat  arbitrary  and  was  selected  primarily  for 
its  ease  of  computation.  Alternative  methods  of  calculating  the  values 
used  for  Q  and  R,  however,  should  cause  minimal  changes  in  the  results. 

Results 

Resonance  Values .  With  the  specification  of  the  necessary  assump¬ 
tions  clearly  stated,  the  first  step  in  obtaining  the  results  was  to 
determine  the  locations  of  the  equilibrium  points  and  the  width  of  the 
resonance  structure  for  each  resonance  for  a  particular  nominal  eccen¬ 
tricity  and  inclination.  Initially  the  structure  was  examined  for  zero 
eccentricity  and  inclination.  Although  the  primary  resonance  band 
displays  structure  for  these  conditions,  examination  of  the  geopotential 
Hamiltonian  for  the  remaining  two  bands  indicates  that  no  structure  will 
result  with  either  zero  eccentricity  or  inclination.  Table  II  lists  the 
results  of  these  conditions  for  the  primary  resonance  (Resonance  1)  but 
not  for  the  remaining  resonances  since  they  will  have  zero  width  and  are 
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of  little  practical  interest.  Positions  of  the  equilibrium  points  are 
given  relative  to  a  nominal  geosynchronous  radius. 

TABLE  II 

Primary  Resonance  Values  —  Zero  Eccentricity  and  Inclination 

Stable  Unstable  Width 

Resonance  1  +  2053.63  m  +  2117.04  m  84313.86  m 

Note  that  the  effect  of  the  nonspherical  nature  of  the  earth's  figure  is 
to  move  the  location  of  both  the  stable  and  unstable  equilibrium  points 
out  by  a  distance  of  approximately  2  kilometers  from  a  nominal  two-body 
geosynchronous  orbit.  This  result  is  in  good  agreement  with  previous 
studies.  More  detailed  results  are  contained  in  Appendix  B. 

To  demonstrate  the  effect  of  a  case  of  non-zero  eccentricity  and 
inclination  of  more  practical  interest  for  use  with  geosynchronous 
communications  satellites,  the  approximate  eccentricity  and  inclination 
for  FltSatCom  V  (International  Designation  1981  073A,  NASA  Catalog 
Number  12635)  were  used.  This  particular  satellite  was  chosen  as  the 
geosynchronous  satellite  having  both  the  largest  eccentricity  and  incli¬ 
nation  listed  in  the  NASA  Satellite  Situation  Report  (Ref  24).  Charac¬ 
teristics  for  the  three  resonance  structures  for  e  -  0.025  and  i  0.11 
radians  are  listed  in  Table  III. 

TABLE  III 

Resonance  Values  —  Non-Zero  Eccentricity  and  Inclination 


Stable 

Unstable 

Width 

Resonance 

1 

+  2012.19  m 

+  2075.03  m 

84192.80 

m 

Resonance 

2 

-  12.87  m 

-  )3.14  m 

246.04 

m 

Resonance 

3 

-  2069.46  m 

-  2069.46  m 

0.28 

m 

Some  obvious  conclusions  may  be  drawn  from  these  results.  First, 
changing  the  eccentricity  and  inclination  does  have  some  effect  on  the 
location  of  the  stable  and  unstable  equilibrium  points  and  the  width  of 
the  resonance  structure,  although  not  very  much.  Second,  with  a  width 
of  over  84  kilometers,  the  primary  resonance  engulfs  the  equilibrium 
points  of  both  the  second  and  third  resonances.  And  finally,  for  any 
reasonable  eccentricity  and  inclination,  the  width  of  the  third  reso¬ 
nance  structure  seems  to  preclude  its  use  for  any  practical  satellite 
placement.  As  a  result  of  this  final  conclusion,  the  third  resonance 
band  was  dropped  from  further  consideration. 

Phase  Portraits .  Once  the  structural  characteristics  of  the  three 
resonance  bands  had  been  determined,  the  next  step  was  to  generate  the 
phase  portraits  of  these  bands.  It  was  decided  that  the  only  phase 
portrait  generated  would  be  that  of  the  primary  resonance.  The  reason¬ 
ing  behind  this  decision  lay  in  the  fact  that  the  general  structure  of 
the  two  resonances  being  considered  was  basically  identical,  with  only 
the  scale  being  different.  Since  the  scale  has  already  been  specified 
in  Table  II,  it  was  felt  to  be  unnecessary  to  generate  more  than  one 
phase  portrait.  Since  the  secondary  objective  of  this  investigation  was 
to  demonstrate  the  feasibility  of  this  technique  as  an  investigative 
tool,  the  decision  to  investigate  other  eccentricities  and  inclinations 
or  to  generate  additional  phase  portraits  was  left  to  future  researchers 
interested  in  more  specific  cases. 

The  results  of  the  generation  of  the  phase  portrait  for  the  primary 
resonance  with  the  values  of  eccentricity  and  inclination  specified  in 
use  for  generating  Table  II  are  illustrated  in  Figure  2.  Due  to 
restrictions  imposed  by  the  scaling  of  the  phase  portrait,  a  conformal 


mapping  was  used  in  the  plotting  of  the  results  rather  than  plotting  the 
results  in  polar  coordinates.  What  appears  to  be  a  straight  line 
slightly  in  from  the  stable  and  unstable  equilibrium  points  is  actually 
the  secondary  resonance's  unstable  equilibrium  contours  drawn  to  the 
same  scale.  The  result  gives  a  feel  for  the  relative  locations  of  the 
two  resonance  structures  and  their  relative  sizes.  The  single  line 
width  is  an  upper  limit  for  the  width  of  this  structure  on  the  scale  of 
the  primary  resonance. 

Librational  Periods.  To  demonstrate  the  more  involved,  and 
presumably  more  precise  method  of  determining  the  period  of  libration 
about  the  stable  equilibrium  points,  the  line  integral  method  was 
applied  to  the  primary  resonance. 

TABLE  IV 

Librational  Periods  —  Primary  Resonance 


Contour 

Period 

Contour 

Period 

Contour 

Period 

1 

670.94  days 

9 

440.81  days 

15 

293.33  days 

440.25  days 

292.77  days 

2 

703.93  days 

10 

398.96  days 

16 

281.00  days 

3 

747.38  days 

398.41  days 

280.44  days 

4 

806.99  days 

11 

368.08  days 

17 

270.13  days 

367.53  days 

269.58  days 

5 

900.46  days 

12 

343.85  days 

18 

260.47  days 

6 

1109.36  days 

343.30  days 

259.91  days 

7 

638.64  days 

13 

324.09  days 

19 

251.79  days 

638.08  days 

323.53  days 

251.24  days 

8 

504.43  days 

14 

307.51  days 

20 

243.94  days 

503.88  days 

306.96  days 

243.38  days 

Table  IV 

lists  the  librational  period  for  each  of 

the  energy 

contours 

shown  in 

Figure  2,  working  from  the 

inside  of  the 

resonance 

structure 

out  (this,  of  course,  excludes  the  unstable  equilibrium  contour  which 
would  have  an  infinite  librational  period).  Use  of  the  line  integral 
method  allows  a  means  for  determining  the  specific  librational  period 
associated  with  a  given  energy  or  displacement  from  the  stable  equilib¬ 
rium  points.  It  should  be  noted  that  decreasing  the  step  size  for  the 
integration  as  a  means  of  verification  produced  minor  changes  (less  than 
one  day)  in  the  resulting  periods. 

Occurrence  of  double  entries  in  Table  IV  are  the  result  of  the 
existence  of  dual  trajectories  for  a  given  energy  level  once  outside  of 
the  stable  equilibrium  region.  The  librational  period  of  each  of  the 
two  contours  of  equal  energy  is  listed,  with  the  outer  of  the  two  bands 
listed  first.  These  periods  may  be  more  appropriately  considered  as  a 
type  of  synodic  period  arising  from  the  relative  drift  of  satellites  in 
different  orbits. 

For  Contours  1  through  6,  which  are  within  the  stable  equilibrium 
region,  the  basic  period  of  libration  is  within  the  range  of  700  to  900 
days  often  seen  cited  in  the  literature  (Ref  1,4,5,6,7,14,16,18).  As  it 
turns  out,  the  period  of  libration  drops  off  quickly  as  a  satellite 
librates  closer  to  the  stable  point  and  approaches  infinity  as  it  nears 
the  boundary  of  the  stable  region.  This  result  is  expected  from  the 
theory  of  unstable  point  analysis. 

Turning  to  the  linearization  method  for  determining  the  period  of 
libration  yields  the  comparable  result  of  667.10  days.  Since  this 
result  is  in  good  agreement  with  the  existing  literature  and  the  results 
of  the  previous  method  and  is  much  less  computationally  intensive  to 
obtain,  this  method  was  used  to  obtain  the  period  of  libration  for  the 
secondary  resonance's  stable  region.  It  should  be  clear  that  it  is 


unnecessary  to  obtain  precise  periods  of  libration  for  either  resonance 
band  for  the  purposes  of  this  study,  that  is,  to  determine  whether  other 
practical  stable  equilibrium  regions  exist  for  the  placement  of  communi¬ 
cations  satellites.  The  result  of  the  determination  of  the  libration 
period  for  the  secondary  resonance  band  is  228,289.46  days  or  approxi¬ 
mately  625  years. 
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Figure  3.  Librational  Period  —  Primary  Resonance 


Combining  the  values  listed  in  Table  IV  with  the  value  obtained 
through  the  linearization  method.  Figure  3  illustrates  how  the  libra¬ 
tional  period  varies  with  changes  in  the  generalized  momentum  S  for  a 


fixed  value  of  the  generalized  coordinate  s  corresponding  to  that  of 
either  of  the  stable  equilibrium  points.  A  minimum  oscillatory  period 
of  667.10  days  occurs  at  the  stable  points.  Moving  radially  away  from 
these  points  results  in  an  rapid  increase  in  the  librational  period  as 
the  unstable  contour  is  approached,  with  the  period  asymptotically 
approaching  infinity.  Once  outside  the  stable  region  indicated  by  the 
two  vertical  lines,  the  period  of  oscillation  quickly  drops  off, 
reflecting  the  relative  drift  of  satellites  in  different  orbits. 

It  should  be  noted  at  this  point  that  the  asstimption  made  in 
Chapter  II  on  the  development  of  this  method  was  verified  numerically. 
That  is,  the  cross  partial  terms  did  indeed  prove  to  be  negligible, 
being  at  least  10  orders  of  magnitude  smaller  than  the  other  term  in 
each  equation. 

Conclusions 

Although  some  of  the  conclusions  to  be  drawn  from  this  study  have 
already  been  noted,  there  remain  several  more  to  be  elucidated.  In 
regard  to  the  first  objective,  it  has  been  demonstrated  that  additional 
stable  equilibrium  points  do  exist  in  the  geosynchronous  regime.  The 
two  being  used  in  the  primary  resonance  band  have  been  well  studied  and 
this  investigation  serves  merely  to  reinforce  previous  conclusions. 

More  significantly,  however,  are  the  results  pertaining  to  the  other  two 
resonance  bands.  Although  the  third  resonance  band  appears  to  be  too 
small  to  be  of  any  practical  use  in  the  cases  considered,  the  second 
resonance  band  appears  to  be  more  than  adequate  in  this  regard. 

The  second  resonance  band  is  not  only  wide  enough  for  placement  of 
communications  satellites,  its  period  of  libration  is  sufficiently  long 
enough  not  to  cause  significant  problems  for  operational  use.  The  only 


tracking  problem  for  ground  stations  for  these  orbits  would  result  from 
the  daily  oscillation  of  the  orbit  due  to  its  inclination. 

The  stable  points  of  the  second  resonance  band  and  their  span  are 
more  than  adequately  separated  spatially  from  those  of  the  primary 
resonance,  and,  it  appears,  can  be  separated  in  longitude  by  adjusting 
the  argument  of  perigee.  If  this  separation  is  indeed  possible,  many 
equilibrium  points  exist,  permitting  placement  of  many  satellites  at  the 
stable  points  within  this  resonance  region  and  allowing  for  a  reduction 
in  station-keeping.  One  problem,  however,  may  prevent  their  practical 
use. 

This  limitation  involves  an  issue  not  directly  addressed  in  this 
investigation.  Since  the  secondary  resonance  band  lies  entirely  within 
the  much  broader  primary  resonance  band,  in  order  to  determine  whether 
or  not  these  new  stable  points  are  of  any  practical  value  would  require 
an  analysis  of  the  combined  effects  of  the  primary  and  secondary  reso¬ 
nance  structures  and  the  amount  of  station-keeping  necessary  to  maintain 
positioning  within  the  secondary  resonance's  stable  region. 

Turning  to  the  secondary  objective,  this  investigation  should  have 
adequately  demonstrated  the  feasibility  of  this  general  technique  as  a 
primary  investigative  tool  for  the  study  of  resonance  structures.  Any¬ 
one  familiar  with  the  analytical  methods  required  for  an  investigation 
of  this  magnitude  will  appreciate  the  relative  simplicity  of  this  numer¬ 
ical  approach  to  the  problem.  It  is  only  hoped  that  this  method  will 
useful  for  further  research  in  this  area. 

Recommendations 

The  primary  avenues  for  further  research  on  this  topic  center  on 
the  limitations  listed  in  the  Conclusions  section  above.  It  is  the 
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author's  belief  that  the  major  thrust  of  future  research  on  this  topic 
should  involve  developing  a  technique  to  analyze  the  stability  of  the 
secondary  resonance's  stable  points  in  conjunction  with  the  primary 
resonance's  effects  in  determining  the  station-keeping  requirements 
necessary  to  maintain  a  satellite  at  one  of  these  points. 

Should  it  prove  feasible  to  maintain  such  orbits,  further  consider¬ 
ation  must  necessarily  be  given  to  the  problem  of  directional  separa¬ 
tion.  It  may  turn  out  that  the  additional  capacity  for  satellite  place¬ 
ment  in  these  orbits  along  with  communications  frequency  separation  will 
provide  for  expanded  use  of  the  geos3mchronous  orbit.  And  finally,  it 
may  be  desirable  to  determine  the  effect  of  the  inclusion  of  additional 
zonal  and  tesseral  harmonic  terms  on  the  conclusions  of  this  and  future 
studies.  Although  such  inclusions  should  not  drastically  alter  the 
overall  structures,  they  may  affect  the  stability  of  these  structures. 

One  further  recommendation  should  also  be  made  in  regard  to  the 
methods  used  for  analytical  techniques.  This  author  found  the  recently 
developed  forms  of  computer  algebra  programs  of  immense  help  in  simpli¬ 
fying  the  tedious  expansions  and  developing  the  partial  derivatives  in 
Chapter  11.  The  use  of  these  programs,  which  are  capable  of  doing 
algebra,  trigonometric  substitutions  and  simplifications,  differentia¬ 
tion,  and  integratio'n,  should  prove  extremely  valuable  in  pushing  for¬ 
ward  the  analytical  methods  without  increasing  the  amount  of  human 
effort  required. 
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Appendix  A 

Computer  Programs  and  Subroutines 


The  following  is  a  listing  of  the  computer  programs  used  in  the 
process  of  this  investigation.  The  four  main  programs  are  listed  first 
followed  by  the  major  subroutines  employed. 


PROGRAM  VALUE 

DOUBLE  PRECISION  J2 , J4 , J22,N0,RAD ,SR,ECC, INC ,X,Y, Z,PI, 

1 STl , UNSTl , RCS , RESMN , RCU, RESMX , DF , ROOTl , R00T2 
INTEGER  RES 

COMMON/S AO/ J2 , J4 , J2 2 , NO , RAD , S R, Y , Z, RES 
ECC  =  2.5D-2 
INC  =  l.lD-1 

Q  ■kirkirk'kii-kleickifkifk-kifkicifk-k-ii-k-k 

C  ***  SAO  III  CONSTANTS  *** 

Q  ■k1c-k1c1e*-k1e1f1rkirk1t1rkirk1t1fkie1e1c-k 

J2  =  1082.637D-6 

J4  =  -1.617999D-6 

J22  =  2.7438636D-6 

NO  *  5.86729371D-2 

RAD  =  6.378140D6 

SR  =  (1 DO/NO )**(2D0/3D0) 

PRINT  '("INOMINAL  RADIUS",F30.26,"  DU"/)',  SR 
X  =  DSQRT(SR) 

Y  =  X*(1D0  -  DSQRTdDC  -  ECC**2)) 

Z  =  X*DSQRT(1D0  -  ECC**2)*(1D0  -  DCOS(INC)) 

PI  =  4D0*DATAN(1D0) 

STl  =  3D0*PI/4D0 
UNSTl  *  PI/4D0 

C  ******************** 

C  ***  MAIN  PROGRAM  *** 

C  ******************** 

DO  1  RES=1,3 

PRINT  '("  RESONANCE  ",I1/)',  RES 
PRINT  '("  STABLE")' 

CALL  RESONCRCS, STl, RESMN) 

PRINT  '("  UNSTABLE")’ 

CALL  RESON(RCU, UNSTl .RESMX) 

DF  «  RESMX  -  RESMN 

CALL  ROOTS ( RCS , S  T1 , DF , ROOTl , ROOT2 , RESMN ) 

PRINT  '("  ROOT  1  -",F30.26)',  ROOTl 
PRINT  '("  ROOT  2  -",F30.26)',  ROOT2 

PRINT  '("  WIDTH  -",F20.5,"  METERS")',  (ROOT2**2  -  ROOTl**2)*RAD 
1  PRINT  '("  DF  -",D34.26/)' ,  DF 
STOP 
END 


The  job  control  for  program  RSPLOT  is  included  due  to  the  machine 


dependent  nature  of  the  plotting  routines  used.  Also  note  that  the 
resonance  of  interest  and  scale  factors  must  be  specified  in  the  FACTORS 
block. 


TSK,T200.  T820730, KELSO, 4289, 10/26/ 82 
ATTACH , CCPLOT , CCPLOTl 038 , ID-LIBRARY, SN=AS D . 

LIBRARY, CCPLOT. 

FTN5. 

LGO. 

ROUTE , TAPE 1 , DC=PU ,T ID-AF , S  T=CS A, FID-TSK . 

*E0R 

PROGRAM  RSPLOT 

DOUBLE  PRECISION  J2 , J4, J22,N0,SR,X,Y, Z,ECC, INC,PI,ST1 ,ST2, 
1 UNSTl , ZE  RO , DF , DFS , DFM , RCS , RESMN , RCU , RESMX , DEL , ROOTl , R00T2 
REAL  XMIN,XMAX,YMIN,YMAX,XSCALE,YSCALE 
INTEGER  RES 

COMMON/SAO/J2,J4, J22,N0,SR,Y,Z,RES 
COMMON/S IZE /XM IN , XMAX , YMIN , YMAX , XS  CALE , YS  CALE 
CALL  PLOTS(0,0,1) 

ECC  =  2.5D-2 
INC  =  l.lD-1 

C  ************************* 

C  ***  SAO  III  CONSTANTS  *** 

C  ************************* 

J2  =  1082.637D-6 

J4  -  -1.617999D-6 

J22  =  2.7438636D-6 

NO  =  5.86729371D-2 

SR  =  (1D0/N0)**(2D0/3D0) 

X  =  DSQRT(SR) 

Y  -  X*(1D0  -  DSQRTdDO  -  ECC**2)) 

Z  -  X*DSQRT(1D0  -  ECC**2)*(1D0  -  DCOS(INC)) 

PI  =  4D0*DATAN(1D0) 

STl  =  3D0*PI/4D0 
ST2  -  7D0*PI/4D0 
UNSTl  -  PI/4D0 
ZERO  -  ODO 

Q  ***************************** 

C  ***  FACTORS  -  RESONANCE  N  *** 

Q  ***************************** 

C  RES  -  N 

C  XMIN  -  ***** 

C  XMAX  -  ***** 

C  DF  -  ***** 

C  DFS  -  ***** 

C  DFM  -  ***** 


RSPLOT  (Continued) 


Q  ******************** 

C  ***  MAIN  PROGRAM  *** 

C  ******************** 

YMIN  -  0.0 

YMAX  *  2.0*SNGL(PI) 

XSCALE  =  (XMAX  -  XMIN)/5.0 
YSCALE  =  (YMAX  -  YMIN)/8.0 
CALL  BOX 

CALL  RES0N(RCS,ST1,RESMN) 

CALL  RES0N(RCU,UNST1,RESMX) 

DEL  =  PI/1D2 

10  IF  (DF  .LT.  DFM)  THEN 

CALL  ROOTS ( RCS , ZE  RO , DF , ROOTl , ROOT2 , RESMN ) 
CALL  DRAW( ROOTl , ZE  RO , DEL ) 

CALL  DRAW( ROOT2 , ZE  RO , DEL ) 

IF  (RESMN  +  DF  .LT.  RESMX)  THEN 

CALL  ROOTS ( RCS , S T1 , DF , ROOTl , ROOT2 .RESMN ) 
CALL  DRAW(ROOT2,STl,DEL) 

CALL  ROOTS ( RCS , S T2 , DF , ROOTl , ROOT2 .RESMN ) 
CALL  DRAW(ROOT2,ST2.DEL) 

CALL  DRAW(ROOT2.ST2,-DEL) 

END  IF 

DF  =  DF  +  DFS 
GO  TO  10 
END  IF 

DF  »=  RESMX  -  RESMN 

CALL  ROOTS ( RCS . ZERO . DF . ROOTl , ROOT2 .RESMN ) 

CALL  DRAW( ROOTl . ZERO . DEL) 

CALL  DRAW(ROOT2. ZERO. DEL) 

CALL  ROOTS ( RCS , S T1 . DF . ROOTl . ROOT2 .RESMN ) 

CALL  DRAW( ROOTl . S  T1 . DEL ) 

CALL  DRAW( ROOTl, STl, -DEL) 

CALL  DRAW(ROOT2,STl,DEL) 

CALL  DRAW(ROOT2,STl,-DEL) 

CALL  ROOTS ( RCS , S T2 , DF , ROOTl , ROOT2 .RESMN ) 

CALL  DRAW( ROOTl, ST2, DEL) 

CALL  DRAW(ROOTl,ST2,-DEL) 

CALL  DRAW(ROOT2,ST2,DEL) 

CALL  DRAW(ROOT2,ST2,-DEL) 

CALL  PLOTE(O) 

STOP 

END 


PROGRAM  TIME 

DOUBLE  PRECISION  J2, J4, J22,NO,SR,ECC, INC,X,Y,Z,PI,ST1,UNST1,RCS, 
IRESMN , RCU , RESMX , DF , DPS , DFM , DEL , ROOTl , ROOT2 
INTEGER  RES 

COMMON/S AO/ J2,J4,J22,NO,SR,Y,Z,RES 
ECC  -  2.5D-2 
INC  -  l.lD-1 

C  ************************* 

C  ***  SAO  III  CONSTANTS  *** 

C  ************************* 

J2  =  1082.637D-6 

J4  *  -1.617999D-6 

J22  =  2.7438636D-6 

NO  =  5.86729371D-2 

SR  =  (1DO/NO)**(2DO/3DO) 

X  =  DSQRT(SR) 

Y  =  X*(1D0  -  DSQRTdDO  -  ECC**2)) 

Z  -  X*DSQRT(1D0  -  ECC**2)*(1D0  -  DCOS(INC)) 

PI  -  4D0*DATAN(1D0) 

STl  =  3D0*PI/4D0 
UNSTl  *  PI/4D0 

C  ***************************** 
c  ***  FACTORS  -  RESONANCE  N  *** 

C  ***************************** 

C  RES  -  N 

C  DF  =  ***** 

C  DFS  »  ***** 

C  DFM  »  ***** 

Q  ******************** 

c  ***  main  program  *** 
c  ******************** 

CALL  RESON(RCS,STl,RESMN) 

CALL  RESON( RCU, UNSTl, RESMX) 

DEL  -  PI/1D3 

10  IF  (DF  .LT.  DFM)  THEN 

CALL  ROOTS ( RCS , S  T1 , DF , ROOTl , ROOT2 , RE  SMN ) 

CALL  PERIOD(ROOT2,STl,DEL) 

IF  (DF  .GT.  RESMX-RESMN)  CALL  PE RIOD( ROOTl , STl , DEL) 

PRINT  '("0")' 

DF  -  DF  +  DFS 
GO  TO  10 
END  IF 
STOP 
END 


PROGRAM  TIMES 

DOUBLE  PRECISION  J2, J4, J22,N0,SR,X,Y,Z,ECC, INC,PI,ST1, 
1RCS,RESMN,XX,YY,DX,DY,DXX,DYY,W,W1,M2,SID,C0N,PER 
INTEGER  RES 

COMMON/SAO/J2 , J4, J22,N0,SR, Y, Z,RES 
ECC  -  2.5D-2 
INC  *  l.lD-1 

C  ************************* 

C  ***  SAO  III  CONSTANTS  *** 

C  ************************* 

J2  =  1082.637D-6 

J4  -  -1.617999D-6 

J22  =  2.7438636D-6 

NO  =  5.86729371D-2 

SR  -  (1 DO/NO )**(2D0/3D0) 

X  -  DSQRT(SR) 

Y  =  X*(1D0  -  DSQRTdDO  -  ECC**2)) 

Z  »  X*DSQRT(1D0  -  ECC**2)*(1D0  -  DCOS(INC)) 

PI  »  4D0*DATAN(1D0) 

SID  -  1.00273790931D0 
CON  -  NO/SID 
STl  *  3D0*PI/4D0 

C  ******************** 

c  ***  main  program  *** 

c  ******************** 

C  RES  =  N 

CALL  RESON(RCS,STl,RESMN) 

DXX  »  lI>-9 

DYY  »  ID-9 

XX  -  RCS  +  DXX 

YY  -  STl  +  DYY 

CALL  FDF(RCS,YY,F,DX,DY) 

W1  -  DSQRT(DABS(DY/DYY)) 

CALL  FDF(XX,ST1,F,DX,DY) 

W2  -  DSQRT(DABS(DX/DXX)) 

W  =  W1  *  W2 
PER  =  COH/W 

PRINT  '("  PERIOD  -",F10.2,"  DAYS")',  PER 

STOP 

END 


Below  is  a  sample  FACTORS  block,  the  one  used  for  this  study. 


Q  ***************************** 

C  ***  FACTORS  -  RESONANCE  1  *** 

Q  ***************************** 

RES  -  1 
XMIN  -  2.57 
XMAX  -  2.58 
DF  -  5D-9 
DFS  -  lD-8 
DFM  -  2D-7 
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The  remainder  of  Appendix  A  is  a  listing  of  major  subroutines 


SUBROUTINE  BOX 

REAL  P(1:7),Q(1:7),XMIN,XMAX,YMIN,YHAX,XSCALE,YSCALE 
COMMON/S IZE/XMIN , XMAX, YMIN , YMAX, XSCALE , YSCALE 
CALL  PLOT(0.0,-0.5,-3) 

CALL  PLOT(0.0,1.5,-3) 

P(l)  »  XMIN 

P(2)  -  XMIN 

P(3)  =  XMAX 

P(4)  =  XMAX 

P(5)  -  XMIN 

P(6)  =  XMIN 

P(7)  =  XSCALE 

Q(l)  =  YMIN 

Q(2)  *  YMAX 

Q(3)  =  YMAX 

Q(4)  =  YMIN 

Q(5)  -  YMIN 

Q(6)  »  YMIN 

Q(7)  -  YSCALE 

CALL  LINE(P,Q,5,1,0,0) 

RETURN 

END 


SUBROUTINE  RESON(S,SY,F) 

DOUBLE  PRECISION  S,SY,F,J2, J4, J22,N0,RAD,SR,Y,Z,DS,DX,DY 
INTEGER  RES 

COMMON/S AO/ J2 , J4 , J2 2 ,N0 ,RAD ,SR, Y , Z , RES 
S  »  2D0 
DS  »  2D0 

30  IF  (DS  .GT.  lD-28)  THEN 
DS  =  DS/2D0 

CALL  FDF(S,SY,F,DX,DY) 

S  »  S  -  DS  *  DSIGNdDO.DX) 

GO  TO  30 
END  IF 

PRINT  '("  RADIUS  »",F30.26,"  DU")’,  S**2 

PRINT  '("  DISTANCE  FROM  NOMINAL  ■",F12.5,"  METERS")', 

1  (S**2  -  SR)*RAD 

PRINT  '("  S  -",F30.26/)' ,  S 
RETURN 
END 
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The  above  listing  of  RESON  was  used  only  in  VALUE.  The  following 
listing  was  used  in  all  other  programs. 


SUBROUTINE  RESON(S,SY,F) 

DOUBLE  PRECISION  S,SY,F,DS,DX,DY 
S  »  2D0 
DS  -  2D0 

30  IF  (DS  .GT.  lD-28)  THEN 
DS  »  DS/2D0 

CALL  FDF(S,SY,F,DX,DY) 

S  -  S  -  DS  *  DSIGNdDO.DX) 

GO  TO  30 
END  IF 
RETURN 
END 


SUBROUTINE  ROOTS (X , Y , DF , ROOTl , ROOT2 , FMIN ) 

DOUBLE  PRECISION  X, Y,DF, ROOTl, ROOT2, FMIN, F,FN, FI, F2,DX,DY 

F  -  FMIN  +  DF 

CALL  FDF(X,Y,FN,DX,DY) 

IF  (FN  ,GT.  F)  THEN 
ROOTl  -  ODO 
ROOT2  -  ODO 
ELSE 

ROOTl  *  2D0 

41  CALL  FDF(R00T1,Y,F1,DX,DY) 

ROOTl  -  ROOTl  -  (F1-F)/DX 

IF  (DABS(F-Fl)  .GT.  lD-28)  GO  TO  41 
ROOT2  -  3D0 

42  CALL  FDF(ROOT2,Y,F2,DX,DY) 

ROOT2  »  ROOT2  -  (F2-F)/DX 

IF  (DABS(F-F2)  .GT.  lD-28)  GO  TO  42 
END  IF 
RETURN 
END 


SUBROUTINE  DRAW(ROOT,Sy,DLT) 

DOUBLE  PRECISION  ROOT,SY,DLT,PI2,DEL,X,Y,XN,YN,F 

REAL  P(l: 2500), Q(l: 2500), XMIN,XMAX,YMIN,YMAX,XSCALE,YSCALE 

INTEGER  N,IOK,ICAS 

COMMON/S IZE /XMIN , XMAX , YMIN . YMAX , XS  CALE , YSCALE 
IF  (ROOT  .EQ.  ODO)  RETURN 
PI2  =  8D0*DATAN(1D0) 

DEL  *  DLT 
N  =  1 

P(N)  *  SNGL(ROOT) 

Q(N)  =  SNGL(SY) 

ICAS  =  0 

CALL  CONTUR( ROOT , S Y , DEL , XN , YN , lOK , ICAS , F ) 

X  =  ROOT 

51  IF  (lOK  .NE.  0)  THEN 

IF  (YN  .GT.  PI2+DEL/2D0  .OR.  YN  .LT.  ODO)  GO  TO  52 

N  =  N  +  1 

P(N)  »  SNGL(XN) 

Q(N)  =  SNGL(YN) 

IF  (DABS(YN-SY)  .LT.  DABS(DEL/2D0)  .AND.  X  .LT.  XN)  GO  TO  52 
X  =  XN 
Y  =  YN 

CALL  CONTUR( X , Y , DEL , XN , YN , lOK , ICAS , F ) 

GO  TO  51 
END  IF 

IF  (N  .EQ.  1)  RETURN 

52  P(N+1)  -  XMIN 
P(N+2)  =  XSCALE 
Q(N+1)  =  YMIN 
Q(N+2)  *  YSCALE 

CALL  LINE(P,Q,N, 1,0,0) 

RETURN 

END 


SUBROUTINE  PERI0D(R00T,SY,DLT) 

DOUBLE  PRECISION  ROOT,SY,DLT,J2, J4, J22,NO,SR,Y,Z,PI2,SID,CONV, 
1DEL,X0,Y0,XN,YN,F,DX,DY,PER 
COMMON/SAO/J2,J4,J22,NO,SR,Y,Z,RES 
INTEGER  RES,IOK,ICAS 
PI2  -  8D0*DATAN(1D0) 

SID  -  1.00273790931D0 

CONV  »  NO/P 12 /SID 

DEL  *  DLT 

XO  =  ROOT 

YO  -  SY 

PER  -  ODO 

ICAS  -  0 

CALL  CONTUR( XO , YO , DEL , XN , YN , lOK , ICAS , F ) 

51  IF  (lOK  .NE.  0)  THEN 

IF  (YN  .GT.  PI2+DEL/2D0)  GO  TO  52 
CALL  FDF(XN,YN,F,DX,DY) 

IF  (ICAS  .EQ.  1)  THEN 

PER  “  PER  +  DABS((YN-YO)/DX) 

ELSE 

PER  »  PER  +  DABS((XN-XO)/DY) 

END  IF 

IF  (DABS(YN-SY)  .LT.  DABS(DEL/2D0)  .AND.  XO  .LT.  XN)  GO  TO  52 
XO  «  XN 
YO  =  YN 

CALL  CONTDR( XO , YO , DEL , XN , YN , lOK , ICAS , F ) 

GO  TO  51 
END  IF 

52  PRINT  '("  PERIOD  -'',F9.2,"  DAYS”)',  PER*CONV 
RETURN 

END 


SUBROUTINE  CONTURCX, Y.DEL.XNEW.YNEW, lOK, ICAS.FF) 

DOUBLE  PRECISION  X,Y, DEL, KNEW, YNEW,FF,TOLXY,TOLF,FDG,F,DFDX,DFDY, 
IDX, DY ,DFPDX,DFPDY , RES , FP .ERRFl  ,ERRF2  ,ERRF .ERRXYl  ,ERRXY2  .ERRXY 
INTEGER  NL,IOK,ICAS,I 
TOLXY  -10-22 
TOLF  -  lD-28 
FDG  -  5D2 
lOK  -  1 
NL  -  50 

CALL  FDF(X,Y,F,DFDX,DFDY) 

IF  (ICAS  .EQ.  0)  FF  =  F 

IF  (DABS(DFDY)  .GT.  DABS(DFDX/FDG) )  THEN 

C  4r4r*4r4r4r4r4r4r4r4r4r4rir'Ar4r4rir4r4r4r 

C  ***  SHALLOW  SLOPE  *** 

C  ********************** 

IF  (ICAS  .EQ.  0)  ICAS  =  -1 
IF  (ICAS  .EQ.  1)  THEN 

DEL  =  -DEL  *  DSIGNdDO.DFDX)  *  DSIGN(1D0,DFDY) 


ICAS  =  -1 


END  IF 

XNEW  =  X  +  DEL/FDG 
YNEW  -  Y 

DY  -  -DFDX*DEL/DFDY/FDG 
DO  61  I  =  1,NL 

YNEW  »  YNEW  +  DY 

CALL  FDF( XNEW, YNEW, FP,DFPDX,DFPDY) 

RES  -  FP  -  FF 
DY  -  -RES/DFPDY 
ERRFl  -  DABS (RES) 

ERRF2  =  DABS(ERRF1/FF) 

ERRF  -  miNl (ERRFl, ERRF2) 

ERRXYl  -  DABS(DY) 

ERRXY2  -  IDO 

IF  (YNEW  .NE.  ODO)  ERRXY2  -  DABS(DY/YNEW) 

ERRXY  =  IWINl (ERRXYl ,ERRXY2) 

IF  (ERRXY  .LT.  TOLXY  .AND.  ERRF  .LT.  TOLF)  GO  TO  63 


61 


CONTUR  (Continued) 


C  ******************* 

C  ***  STEEP  SLOPE  *** 

Q  ******************* 


ELSE 

IF  (ICAS  .EQ.  0)  ICAS  =  1 
IF  (ICAS  .EQ.  -1)  THEN 

DEL  =  -DEL  *  DSIGNdDO.DFDX)  *  DS IGN ( 1  DO , DFDY  ) 

ICAS  =  1 
END  IF 
XNEW  =  X 
YNEW  =  Y  +  DEL 
DX  *  -DFDY*DEL/DFDX 
DO  62  I  -  1,NL 

XNEW  -  XNEW  +  DX 

CALL  FDF(XNEW,YNEW,FP,DFPDX,DFPDY) 

RES  =  FP  -  FF 
DX  =  -RES/DFPDX 
ERRFl  «  DABS (RES) 

ERRF2  =  DABS(ERRF1/FF) 

ERRF  =  miNl  (ERRFl,  ERRF2) 

ERRXYl  =  DABS(DX) 

ERRXY2  «  DABS(DX/XNEW) 

ERRXY  =  DMINl (ERRXYl, ERRXY2) 

62  IF  (ERRXY  .LT.  TOLXY  .AND.  ERRF  .LT.  TOLF)  GO  TO  63 


END  IF 
lOK  =  0 
63  RETURN 
END 


nnn  000  000 


The  COMMON  and  DOUBLE  PRECISION  blocks  of  PDF  should  contain  the 


variable  RAD  when  run  with  program  VALUE. 


SU  BROUT INE  FDF ( X , XX , F , DFDX , DFDXX ) 

DOUBLE  PRECISION  X, XX, F, DFDX, DFDXX, J2,J4,J22,N0, RAD, SR, Y,Z,Q,R,S, 
lA , B , C , FO , F 1 , F2 , F4 , R22 , DFDXO , DFDXl , DFDX2 , DFDX4 , DFDX22 
INTEGER  RES 

COMMON/S AO/ J2 , J4 , J2 2 , NO , RAD , S R , Y , Z , RES 
Q  »  Y  -  DBLE(RES-1)*DSQRT(SR) 

R  =  Z  +  DBLE(RES-1)*DSQRT(SR) 

S  =  X 

GO  TO  (71,72,73)  RES 

C  ********************* 

C  ***  1ST  RESONANCE  *** 

C  ********************* 

71  A  =  S/(S  -  Q) 

B  =  R/(S  -  Q) 

C  =  Q/S 

**************** 

***  J22  TERM  *** 

**************** 

R22  -  J22/(X**6)*(4D0  -  4D0*B  +  B**2)*(288D0  -  144D1*C 
1+  1656D0*C**2  -  1216D0*C**3  +  654D0*C**4  -  21D1*C**5  +  35D0*C**6) 
2*DSIN(2D0*XX)/384D0 
*************************** 

***  momentum  derivative  *** 

*************************** 

DFDXl  =  J2/(X**7)*(A**3*(-  6D0  +  18D0*B  -  9D0*B**2)  +  A**4*(-  6D0 
1+  24D0*B  -  15D0*B**2))/4D0 

DFDX2  »  J2**2/(X**11)*(A**5*(12D1  -  24D1*B  -  18D1*B**2 
1+  3D2*B**3  -  75*B**4)  +  A**6*(-  72D0  +  864D0*B  -  2556D0*B**2 
2+  2208D0*B**3  -  567D0*B**4)  +  A**7*(-  648D0  +  3456D0*B 
3-  6588D0*B**2  +  5148D0*B**3  -  1395D0*B**4)  +  A**8*{-  84D1 
4+  384D1*B  -  594D1*B**2  +  42D2*B**3  -  1155*B**4))/128D0 
DFDX4  *  3D0*J4/(X**11)*((A**5  -  A**7)*(-  12D1  +  12D2*B 

1-  27D2*B**2  +  21D2*B**3  -  525D0*B**4)  +  A**6*(-  12D1 
2+  144D1*B  -  378D1*B**2  +  336D1*B**3  -  945D0*B**4)  +  A**8*(28D1 
3-  32D2*B  +  81D2*B**2  -  7D3*B**3  +  1925D0*B**4) )/128D0 

DFDX22  =  3D0*J22/(X**7)*((4D0  -  4D0*B  +  B**2)*(-  288D0 
1+  168D1*C  -  2208D0*C**2  +  1824D0*C**3  -  109Dl*C**4  +  385D0*C**5 

2-  7D1*C**6)/48D0  +  A*(4D0*B  -  2D0*B**2)*(288D0  -  144D1*C 
3+  1656D0*C**2  -  1216D0*C**3  +  654D0*C**4  -  21D1*C**5  +  35D0*C**6) 
4/ 288D0 )*DS IN( 2D0*XX )/4D0 

*************************** 

***  POSITION  DERIVATIVE  *** 

*************************** 

DFDXX  -  J22/X**6*(4D0  -  4D0*B  +  B**2)*(288D0 
1-  144D1*C  +  1656D0*C**2  -  1216D0*C**3  +  654D0*C**4  -  21D1*C**5 
2+  35DO*C**6)*DCOS(2DO*XX)/192DO 
GO  TO  74 


oooooo  ooo 


FDF  (Continued) 


C  ********************* 

C  ***  2ND  RESONANCE  *** 

C  ********************* 

72  A  =  -S/Q 

B  =  (S  -  R)/Q 
C  =  (Q  +  S)/S 
**************** 

***  J22  TERM  *** 

**************** 

R22  =  J22/(X**6)*(3D0/128D0*(2D0*B  -  B**2)*(288D0*C  +  304D0*C**2 
1+  68D1*C**3  -  158D1*C**4  +  846D0*C**5  -  141D0*C**6) )*DSIN(2D0*XX) 
*************************** 

***  MOMENTUM  DERIVATIVE  *** 

*************************** 

DFDXl  =  J2/(X**7)*(A**3*(-  6D0  +  18D0*B  -  9D0*B**2)  +  A**4*(6D0 
1-  6D0*B))/4D0 

DFDX2  =  J2**2/(X**11)*(A**5*(12D1  -  24D1*B  -  18D1*B**2 
1+  3D2*B**3  -  75D0*B**4)  +  A**6*(-  24D1  +  108D1*B  -  2124D0*B**2 
2+  1668D0*B**3  -  432D0*B**4)  +  A**7*(-  72D0  +  288D0*B  -  684D0*B**2 
3+  828D0*B**3  -  315D0*B**4)  +  A**8*(48D1  -  132D1*B  +  126D1*B**2 
4-  42D1*B**3))/128D0 

DFDX4  =  3D0*J4/(X**11)*((A**5  -  A**7)*(-  12D1  +  12D2*B 

1-  27D2*B**2  +  21D2*B**3  -  525D0*B**4)/1 6D0  +  (3D0*A**6 

2-  5D0*A**8)*(6D1  -  27D1*B  +  315D0*B**2  -  105D0*B**3)/12D0)/8D0 
DFDX22  *  3D0*J22/(X**7)*((2D0*B  -  B**2)*(144D0  -  704D0*C 

1-196D0*C**2  -  622D1*C**3  +  10015D0*C**4  -  5076D0*C**5 
2+  846D0*C**6)/64D0  -  A*(2D0  -  2D0*B)*(288D0*C  +  304D0*C**2 
3+  68D1*C**3  -  158D1*C**4  +  846D0*C**5  -  141D0*C**6)/128D0) 
4*DSIN(2D0*XX) 

*******  *******************4r 

***  POSITION  DERIVATIVE  *** 

*************************** 

DFDXX  =  J22/(X**6)*(3D0/64D0*(2D0*B  -  B**2)*(288D0*C  +  304D0*C**2 
1+  68D1*C**3  -  158D1*C**4  +  846D0*C**5  -  141D0*C**6) )*DCOS(2DO*XX) 
GO  TO  74 


oooooo  uuu 


FDF  (Continued) 


C  ********************* 

C  ***  3RD  RESONANCE  *** 

C  ********************* 

73  A  *  -S/(S  +  Q) 

B  *  (2D0*S  -  R)/(S  +  Q) 

C  =  (2D0*S  +  Q)/S 
**************** 

***  J22  TERM  *** 

**************** 

R22  =  J22/(X**6)*B**2*(4D1*C**2  +  16D0*C**3  -  74D0*C**4 
1+  42D0*C**5  -  7D0*C**6)*DSIN(2D0*XX)/32D1 
*************************** 

***  MOMENTUM  DERIVATIVE  *** 

*************************** 

DFDXl  =  J2/(X**7)*(A**3*(-  6D0  +  18D0*B  -  9D0*B**2)  +  A**4*(18D0 

1-  36D0*B  +  15D0*B**2))/4D0 

DFDX2  =  J2**2/(X**11)*(A**5*(12D1  -  24D1*B  -  18D1*B**2 
1+  3D2*B**3  -  75D0*B**4)  +  A**6*(-  408D0  +  1296D0*B  -  1692D0*B**2 
2+  1128D0*B**3  -  297D0*B**4)  +  A**7*(504D0  -  288D1*B  +  522D1*B**2 
3-  3492D0*B**3  +  765D0*B**4)  +  A**8*(18D2  -  648D1*B 
4+  846D1*B**2  -  504Dl*B**3  +  1155D0*B**4) )/128D0 
DFDX4  =  3D0*J4/(X**11)*((A**5  -  A**7)*(-  12D1  +  12D2*B 
1-  27D2*B**2  +  21D2*B**3  -  525D0*B**4)  +  A**6*(6D2  -  36D2*B 
2+  63D2*B**2  -  42D2*B**3  945D0*B**4)  +  A**8*(-  108D1 

3+  68D2*B  -  123D2*B**2  +  84D2*B**3  -  1925D0*B**4) )/128D0 
DFDX22  =  3D0*J22/(X**7)*(B**2*(8D1*C  -  112D0*C**2  -  368D0*C**3 
1+  58D1*C**4  -  273D0*C**5  +  42D0*C**6)  -  A*(2D0*B  -  B**2) 
2*(4D1*C**2  +  16D0*C**3  -  74D0*C**4  +  42D0*C**5  -  7D0*C**6)) 
3*DSIN(2D0*XX)/48D1 
*************************** 

***  POSITION  DERIVATIVE  *** 

*************************** 

DFDXX  =  J22/(X**6)*B**2*(4D1*C**2  +  16D0*C**3  -  74D0*C**4 
1+  42D0*C**5  -  7DO*C**6)*DCOS(2DO*XX)/16D1 

C  ******************* 

C  ***  HAMILTONIAN  *** 

Q  ******************* 

74  FO  =  1D0/(X**2)/2D0 

FI  »  J2/(X**6)*A**3*(2D0  -  6D0*B  +  3D0*B**2)/4D0 

F2  »  J2**2/(X**10)*(3D0/32D0*A**5*(-  8D0  +  16D0*B  +  12D0*B**2 

1-  2D1*B**3  +  5D0*B**4)  +  3D0/8D0*A**6*(4D0  -  24D0*B  +  48D0*B**2 

2-  36D0*B**3  +  9D0*B**4)  -  15D0/32D0*A**7*(-  8D0  +  32D0*B 

3-  44D0*B**2  +  28D0*B**3  -  7D0*B**4) )/4d0 

F4  -  J4/(X**10)*(9D0*A**5  -  15D0*A**7)*(8D0  -  8D1*B 
1+  18D1*B**2  -  14D1*B**3  ♦  35D0*B**4)/128D0 
F  -  FO  +  FI  +  F2  +  F4  +  R22  +  N0*X 
DFDXO  -  -  1D0/X**3 

DFDX  -  DFDXO  +  DFDXl  +  DFDX2  +  DFDX4  +  DFDX22  +  NO 

RETURN 

END 


Appendix  B 
Computer  Output 


Appendix  B  contains  a  summary  of  the  output  generated  for  the 
program  VALUE.  All  other  relevant  computer  output  is  summarized  within 
the  text  of  the  thesis.  The  first  set  of  values,  for  Resonance  1  only, 
was  computed  for  zero  eccentricity  and  inclination.  All  other 
calculations  assume  e  =  0.25  and  i  =  0.11  radians. 


NOMINAL  RADIUS  6.62279705974354939872109534  DU 

RESONANCE  1 

STABLE 

RADIUS  =  6.62311903956047332338676051  DU 

DISTANCE  FROM  NOMINAL  »  2053.63235  METERS 

S  *  2.57354211925130017355733216 

UNSTABLE 

RADIUS  >>  6.62312898069230474566012610  DU 

DISTANCE  FROM  NOMINAL  -  2117.03828  METERS 

S  *  2.57354405066093724823231373 

ROOT  1  »  2.57225535879742894051568446 

ROOT  2  *  2.57482973825684011213949358 

WIDTH  »  84513.86393  METERS 

DF  =  .56666332874386549500224905D-07 


RESONANCE  1 
STABLE 

RADIUS  -  6.62311254138901482663038207  DU 

DISTANCE  FROM  NOMINAL  -  2012.18610  METERS 

S  »  2.57354085675534104714050865 

UNSTABLE 

RADIUS  -  6.62312239460014021257026258  DU 

DISTANCE  FROM  NOMINAL  «  2075.03126  METERS 

S  -  2.57354277108427716384806879 

ROOT  1  -  2.57225898236003338751447668 

ROOT  2  -  2.57482358318862921027817112 

WIDTH  -  84192.80152  METERS 

DF  -  .56236552674965968019912343D-07 
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RESONANCE  2 


STABLE 

RADIUS  “  6.62279504161728455154854446  DU 

DISTANCE  FROM  NOMINAL  =  -12.87189  METERS 

S  =  2.57347917062044326769576973 

UNSTABLE 

RADIUS  6.62279499906509052836411341  DU 
DISTANCE  FROM  NOMINAL  =  -13.14330  METERS 

S  =  2.57347916235299844488103640 

ROOT  1  =  2.57347542322313551926585699 

ROOT  2  =  2.57348291802502672216083357 

WIDTH  =  246.03951  METERS 

DF  =  .48024587235371029523549358D-12 

RESONANCE  3 

STABLE 

RADIUS  »  6.62247259764080133997643103  DU 

DISTANCE  FROM  NOMINAL  =  -2069.46472  METERS 
S  =  2.57341652237658591978112241 

UNSTABLE 

RADIUS  -  6.62247259764079875539562401  DU 

DISTANCE  FROM  NOMINAL  =  -2069.46472  METERS 
S  =  2.57341652237658541761196624 

ROOT  1  =  2.57341652195460740223585979 

ROOT  2  =  2.57341652279856443771039969 

WIDTH  -  .02770  METERS 

DF  =  .60909469144020180570086397D-20 


vita 
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